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Abstract 

In the context of the recently developed "equation-free" approach to the computer- 
assisted analysis of complex systems, we illustrate the computation of coarsely self- 
similar solutions. Dynamic renormalization and fixed point algorithms for the macro- 
scopic density dynamics are applied to the results of short bursts of appropriately 
initialized molecular dynamics in a simple diffusion simulation. The approach holds 
promise for locating coarse self-similar solutions and the corresponding exponents in a 
variety of multiscale computational contexts. 
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1. Introduction 

In contemporary scientific and engineering modeling we are often faced with sit- 
uations where the best system models are available at a fine scale (e.g. atomistic, 
individual-based); yet we want to predict the behavior of the system at a much more 
coarse-grained, macroscopic level. In traditional modeling successful closures often al- 
low us to write models directly at the coarse-grained, macroscopic level at which we 
want to model the behavior; typical examples include chemical kinetics closures in 
terms of reactant concentrations for reactor modeling, or Newton's law of viscosity in 
the Navier-Stokes equations. Often, however, such accurate closures are not available, 
and the immense range of active scales (in space as well as time) precludes the practical 
prediction of macroscopic behavior through direct atomistic simulation. 

The recently developed "equation free" approach to coarse-grained computer-assisted 
analysis of complex systems attempts to bridge this enormous scale gap when macro- 
scopic equations conceptually exist, but are not available in closed form jll2l3l4l5l6l7l8j . 
The approach constitutes a bridge between traditional continuum numerical analysis 
and microscopic simulation. The main idea is to consider the microscopic simulation as 
a computational experiment which can be initialized and run at will. If a coarse-grained, 
macroscopic equation were available, traditional numerical tasks would involve repeated 
evaluations of the equation, and of its functional or parametric derivatives, at various 
values of the macroscopic state. The idea is to substitute these function evaluations 
with short bursts of appropriately initialized microscopic computational experimenta- 
tion, from the results of which the requisite numerical quantities are estimated. The 
quasi-continuum method of Phillips, Ortiz and coworkers ,9|, as well as the optimal 
predictor approach of Chorin and coworkers ^10^ embody many of these features; see [7] 
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for a discussion. The approach can be successfully combined with matrix-free iterative 
linear algebra to allow us to solve linear and nonlinear macroscopic equations, perform 
"coarse projective integration" as well as additional tasks like controller design and op- 
timization without ever writing down the macroscopic equations in closed form. This 
is a system identification based, "closure on demand" approach. 

Over the last few years we have demonstrated how to use this approach for the accel- 
erated simulation and bifurcation analysis of the coarse-grained, expected behavior of 
kinetic Monte Carlo, Brownian Dynamics, Molecular Dynamics and Lattice-Boltzmann 
microscopic simulation codes, 01211211^131011313^2 as well as the solution of effective 
medium equations jl2lll3j for media with spatially varying properties. The approach 
gives rise to two-level (conceptually possibly multi-level) codes. The "inner code" is 
the best microscopic simulation of the phenomenon at our disposal. The outer code, 
the "wrapper", is typically templated on traditional continuum numerical analysis, and 
constitutes a protocol for the design of microscopic simulations and for processing their 
results towards a macroscopic modeling goal. Accelerated simulation, the location of 
coarse fixed points and their continuation and bifurcation analysis are typical such 
goals. 

In many cases of interest, the macroscopic dynamics we explore do not involve 
stationary solutions, but rather traveling or self-similar solutions. In this paper we 
will show how to use the basic coarse timestepper methodology to construct dynamic 
renormalization algorithms for the location of coarse self-similar solutions and the 
corresponding similarity exponents by acting on a microscopic code directly. Our illus- 
trative example - molecular diffusion in a thin two dimensional domain - is extremely 
simple, yet it does illustrate the basic ingredients of the approach. 
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The paper is organized as follows: In the next Section we will briefly describe 
template-based dynamic renormalization when macroscopic equations are explicitly 
available. In Section 3 we summarize the main components of the coarse timestepper, 
and, in particular, its modification for computation of self-similar solutions as fixed 
points of the appropriate discrete time map. Section 4 presents our illustrative ex- 
amples for both coarse dynamic and coarse fixed point computations in a dynamic 
renormalization context. We conclude with a brief summary and outline of the scope 
of the method and some of the challenges we expect to arise in its wider application. 

2. Dynamic renormalization in a continuum context 

In problems with translational invariance, one often encounters traveling wave so- 
lutions - constant shape waves that move in space at constant speed. It is convenient 
(mathematically, computationally and practically) to study these solutions, and the 
transient approach to them, in a co-traveling frame. In this frame the traveling so- 
lution appears stationary, and it is much easier to study the transient approach to 
it and its stability unencumbered by its constant motion. Good techniques exist for 
computationally locating the translationally invariant solution along with its traveling 
speed as a nonlinear eigenvalue problem. During transient simulation, however, the 
solutions both travel and approach their ultimate, translationally invariant shape; the 
right speed for "moving along" with such a transient solution may change from moment 
to moment, and the best way to choose it is not transparent. In a recent paper Rowley 
and Marsden described a template based approach that allows one to systematically 
recover such an appropriate instantaneous speed; as the transient solution asymptoti- 
cally approaches the ultimate, translationally invariant, traveling wave, the speed from 
the template-based algorithm asymptotically approaches the correct traveling speed. 
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In a sequence of papers |15lll6lll7j we have shown how to adapt this approach from 
the computation of translationally invariant solutions to the case of scale invariant 
solutions - that is, solutions of dynamic equations that evolve across scales. Self- 
similar solutions |181ll9j are an important such class of scale-invariant solutions; in the 
same sense that it is practical to observe a traveling solution in a co-traveling frame, 
self-similar solutions are convenient to observe in a co-exploding or co-collapsing frame. 
Consider the rather general form of partial differential equation 



for ^ > 0, > 0, and a and b are constants. We assume that there exists a family of 
self-similar solutions. 



where s = \t — t*\ (for the case of finite time blow-up at t = t*) or s = t otherwise, and 
U satisfies the ODE, 




(1) 



where the operator D satisfies the scaling relation 



D^{Bf{x/A)) = A'^B'Dyifiy)), where y = x/A 



(2) 




(3) 



apU - ayUy = Dy{U) 



(4) 



where y = x/{s") with a = sgn{t — t*) and 



P - 1 = aa + b(3. 



(5) 
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We are interested in locating the self-similar shape of the solution as well as both 
similarity exponents a and /?: one more equation is needed for this latter task. Starting 
with the general scaling 



where A, B and r are unknown functions, and setting Ts(s) = aA°'B^ ^, the PDE 
becomes 



This is the "co-exploding" or "co-collapsing" equation which, for self-similar problems, 
is analogous to the "co-traveling" equation for translationally invariant ones. Two 
additional constraints, frequently called pinning conditions, are needed to find both 
A{t) and B{t). There is some latitude in the selection of these conditions - differ- 
ent conditions correspond to different ways of "moving along with the solution across 
scales" before Eq. (O reaches steady state, although all appropriate pinnings give the 
same self-similar shape and exponents 20.,. In the spirit of Rowley and Marsden, we 
proposed in ^S] that such conditions can be constructed by imposing relationships 
between the solution and (essentially arbitrary) template functions. 

Since there are two degrees of freedom ("amplitude" and "width") we must impose 
two pinning conditions. Once these have been imposed, we can use Eq. Q to solve 
for A[t) and B{t): it is actually possible to eliminate the At/ A and B^/B terms in 
Eq. ^ to end up with a "co-exploding" PDE that we called in ^3] MN-dynamics. 
When (and if) the solution of this PDE approaches an asymptotic steady state, we can 




(6) 
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compare coefficients in Eqs. (01) and (O to find: 



(8) 



/3 lim^^oo-Br/^ 



Eqs. © and Q allow us to obtain the scaling exponents a and /3. A more detailed 
discussion can be found in 20 ; the approach can be used to locate both types of self- 
similar solutions JH]) and indeed in ^S] it was used to locate both the Barenblatt and 
the Graveleau solutions of the porous medium equation. 

In view of our illustrative example using molecular dynamics below, we briefly study 
the very simple case of the 1-dimensional diffusion equation, 



The operator in this case is Dx{u) = Uxx- From Eq. 0, 6 = 1 and a = —2. Eq. 
then yields a = 1/2. For self-similar solutions of the second kind, a and (3 cannot 
be computed a priori, and will be found as part of the solution process. For the two 
pinning conditions, we choose here to let w satisfy the following relations to two given 
template functions, Ti{y) and T2{y) respectively, 



These can be thought of as controlling the "width" and the amplitude of the solution. 
By multiplying Eq. ((7| with each T{y), integrating and using the above two equations, 
we can eliminate A and B. A discussion of important technical conditions, having 



(9) 




= 



(10) 




(11) 
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to do with the existence of these integrals over infinite domains, and the appropriate 
solution spaces, can be found in j34ll35j . (An alternative approach is to require that the 
difference between the solution of Eq. Q and a template is minimized, as discussed 

in [laini.) 

For our simple diffusion example, we set Ti{y) to be 

1 for \y\ < 1/2 
Ti{y) = { (12) 

-1 for \y\ > 1/2. 

From Eqs. 0, (fTU]) . and (fT^ we get an equation for Ar/A. Since the subspace of 
solutions symmetric around x = is invariant, restricting our search to symmetric 
solutions we find that 

Ar _ 2Wy{0.5,T) 



A w{0.5,t) 
Similarly, if we set T2{y) = 6{y) we get 



(13) 



Br 

B 



Wyy{f),T). (14) 



Substituting these in Eq. Q we get 



, , 2'u;„(0.5,r) ^ ^ 

Wr+Wyy{Q,T)W + ^^^^—^yWy=Wyy. (15) 



For an initial condition consistent with the constraints and the symmetry we take 



1 for \y\ < 1 

w{y,^) = { (16) 
for \y\ > 1. 
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so that both Eqs. H1U|) and are satisfied. Finally, we numerically integrate Eq. 
((T3|) to steady state, i.e. r — oo, and evaluate the exponent, f3 from 



Another choice for T2{y) would be T2{y) = 1. Physically it corresponds to conser- 
vation of mass, and the same procedure leads to 



Without any further integration, it is clear that /3 = —a = —1/2. 

The above methodology evolves the differential equation in a dynamically rescaled 
time and space frame. Figure^ shows the solution of the diffusion equation (in the 
central portion of a large domain) for the initial condition above; as time progresses we 
know that the solution decays, but it also asymptotically approaches a (self-similar ly 
decaying) Gaussian. (For clarity, the solution is shown over several small blocks of 
time, alhough the equation was also integrated over the "spaces" in the figure.) Figure 

shows the same evolution in a rescaled frame using Eq. (|15|) : the decay has now 
been removed through the template-based rescaling, and one only sees the transient 
approach to the self-similar shape (the Gaussian) consistent with the template-based 
pinning conditions. Figure^ shows the same results in terms of cumulative material 
density and not the density itself; as we will discuss below, it is numerically more 
convenient in particle based simulation codes to work with the cumulative distribution 
function rather than the particle distribution function itself. It is interesting to consider 
the case in which we have a so-called "legacy code" - a code that evolves the original 



aw{0.5, +oo)wyy{0, +00) 



(17) 



2wy{0.5,+oo) 




(18) 
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equation and which we can run but cannot modify. The so-called numerical analysis of 
legacy codes allows us to transform a direct legacy dynamic simulator, by '"wrapping" 
a computational superstructure around it, into a code capable of performing a different 
set of tasks, for which the legacy simulator was not designed. In the dynamic renor- 
malization case, we will compute self-similar solutions by evolving in physical variables 
and rescaling the results, as opposed to first obtaining and then evolving the rescaled 
equation. Our approach is a discrete-time approach (see also [T^21j ^: pioneering work 
on dynamic renormalization, especially in the context of the Nonlinear Schrodinger 
Equation can be found in [^OHll^HOnil^lTTll^l^rMrTT] . 

3. Timestepping for coarse dynamic renormalization 

As we discussed above, one can use a direct simulator of the original equation to 
converge computationally to a (stable) self-similar solution as follows. Starting with 
an initial profile, one evolves forward for a finite time; one then uses the template 
conditions to rescale the space variable for the final profile. The rescaled profile is 
then given as an initial condition to the original equation, which is again evolved 
over finite time; the space variable is rescaled, and the procedure is iterated until 
the shape converges to a member of the family of self-similar shapes. The idea here is 
that rescaling, and then evolving the rescaled equation for a finite time, commutes with 
evolving in physical space and then rescaling the result (see Figure|21). Also, while stable 
self-similar solutions can be found by such dynamic rescaling and forward integration, 
they can also be found - and so can unstable self-similar solutions - through fixed point 
algorithms (like Newton- Raphson) . Indeed, if we call ^Tiw{y)) the result of integrating 
the rescaled equation with initial condition w{y) over time T, the self-similar solution 
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satisfies 

w{y) - <pTiw{y)) = 0. (19) 

We have shown in the past how matrix-free iterative linear algebra techniques can 
be used to converge to solutions of such problems even when the only available tool 
is a subroutine that numerically computes ^Tiw{y)). The original inspiration for this 
work was the so-called Recursive Projection Method (RPM) of Shroff and Keller |33j . 
who used this subroutine (the timestepper) and a computational superstructure (the 
RPM wrapper) to construct a fixed point algorithm. Under appropriate conditions 
this algorithm accelerates the computation of stable fixed points and also stabilizes 
the computation of dynamically unstable ones. In effect, timestepping combined with 
matrix-free techniques "fools" a dynamic simulator into becoming a fixed point solver. 
It is important to note in this entire exposition that we have essentially ignored here 
the role of the boundary conditions for the original and the rescaled equation, assum- 
ing we are solving in sufficiently large domains and for sufficiently short times; it is 
conceivable that weighted Sobolev spaces must be used in order to avoid spurious nu- 
merical oscillations This is an important issue which must be studied carefully; 
yet, as we will see, for our simple problem this does not create major computational 
difficulties. 

We now return to the premise of our introduction: we have a microscopic code 
(e.g. molecular dynamics, evolving a distribution of molecules); yet we believe that 
the coarse-grained, macroscopic behavior of the statistics of the simulation satisfies a 
macroscopic equation that possesses self-similar solutions. We will find these solutions 
through what we will call the "coarse timestepper" , which we have extensively discussed 
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in ^miSllSElinilTllHl^^, and which - for the case of dynamic renormaUzation - is 
ihustrated in Figure 01 This "coarse dynamicahy renormalized timestepper" consists 
of the fohowing elements: 

1. Choose the statistics of interest for describing the coarse-grained behavior of the 
system and an appropriate representation for them. In this case the concentration 
profile in one space dimension is the appropriate macroscopic observable; it is 
the zeroth moment of the distribution of molecules over velocities (and over the 
second, "thin" dimension). It is more convenient to use the particle instantaneous 
Cumulative Distribution Function; assuming it is smooth enough, one can use a 
low-dimensional description of it based on the first few of an appropriate sequence 
of orthogonal polynomials j36ll37j . We will call this the macroscopic description, 
u. These choices determine a restriction operator, A4, from the microscopic-level 
description, U (particle coordinates) to the macroscopic description: u = Ai\J. 

2. Choose an appropriate lifting operator, fi, from the macroscopic description, u, to 
the microscopic description, U. In our case we make random particle assignments 
consistent with the macroscopic CDF. fi should have the property that Aifi is 
the identity (A^/i = I)- In other words, lifting from the macroscopic to the 
microscopic and then restricting (projecting) down again should have no effect 
(except round-off). 

3. From an initial value at the microscopic level, \J{tQ), run the microscopic simula- 
tor (the fine timestepper) for a (relatively short) macroscopic reporting horizon 
T generating the values U(r). We may have to repeat this for several micro- 
scopic initial conditions \Ji{tQ), consistent with the same macroscopic one u(to), 
for variance reduction purposes. 



4. Obtain the restriction u(T) = A4U(T) (the average restriction, in the case of 
many copies). 

5. Rescale u(r) to obtain ur(T) (using the template conditions). 

6. Lift UR,(r) to get a new consistent microscopic U = ^ur,(T) and use it as a new 
starting value for repeating steps 3 to 6. 

Since the diffusion equation has a stable self-similar solution, one can simply repeat 
the above procedure and observe the approach of the statistics of the molecular descrip- 
tion to the Gaussian; the repeated dynamic coarse rescaling helps avoid the continuous 
decay of the direct simulation towards zero. Alternatively, as we will show, one can 
use coarse fixed point algorithms (such as Newton-Raphson) to converge iteratively to 
the fixed point of the coarse rescaled timestepper (as opposed to repeatedly integrating 
and rescaling). Finally, while this is a problem where the exponents are known a priori 
through scaling, the formulas presented in Section 2 and the rescaling history upon 
convergence to the self-similar solution helps us estimate the limiting values of Ar/A 
and B-t/B. This will then help recover, through Eqs. © and © the self-similarity 
exponents for type-2 self-similar solutions. 

4. The computational experiment 

In this Section we outline the data collection procedure from the molecular dynamics 
simulation. A standard Lennard-Jones potential was used, i.e. 




(20) 



with cutoff distance 2.5fT. In what follows all results will be given in reduced units, 
i.e., length in units of cr, time in units of (mcr^/e)^/^. The simulation was performed 
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in a quasi-ld box with size 400 x 10, and periodic boundary conditions in both x and 
y directions. We use T = 1.0, p = 0.5. The domain for such a simulation time was 
"large enough" to appear infinite over the time of our simulation. 

First, the system is allowed to evolve to thermal equilibrium (as evidenced by 
stationarity of the average kinetic energy of the particles). After some simulation time 
we record our first equilibrium configuration, and "reset the clock" to t = 0. Subsequent 
equilibrium configurations are recorded at later times. To model the diffusion process, 
a fraction of the particles were colored according to prescribed distributions (consistent 
with given density profiles), and their positions in the later configurations were tracked. 

We are going to work with a single- variable cumulative distribution function (CDF) 
[3f)] . Suppose that particles are colored, and the i-th particle is in position Xi. Given 
all the colored particle positions, we can immediately obtain their empirical CDF by 
sorting them, that is, relabeling them so that Xj > Xi^i, < i < N , and plotting Xi vs 
Pi = (i — 0.5)/N). We use the CDF rather than the density function because it is trivial 
to generate the empirical CDF while it is computationally difficult and error-prone to 
compute the empirical density function which is the (ill-defined) derivative of the CDF. 

Using only the x coordinate of particle positions is justified by our quasi-ld simu- 
lation box. In fact, it is easier to work with the inverse CDF, or ICDF, Q{p)- It gives 
the the x coordinate of a given particle position, i, and it can be read off directly as 
Q{p{i))- Such continuous ICDFs can then be microscopically approximated by color- 
ing the particle with the nearest x-coordinate value in the simulation box. Using the 
data (positions of colored particles) from later time configurations, we can estimate 
the dynamics of the CDF (or ICDF) evolution. 

In the first numerical experiments, we ran the systems for short bursts of time, then 
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applied templates to rescale the ICDF. Because the number of colored particles we used 
is constant, mass is automatically conserved, so that only one additional template, Ti, 
is needed for dynamic rescaling. For this second template, we computed the slope of 
the center 20% of the colored particles by linear least squares, and normalized it to 
a fixed value. This implied a simple rescaling of the x coordinates of the computed 
colored particle positions. These rescaled x values were then used to "lift" the rescaled 
ICDF: we selected the closest particles to these coordinates in the full set, and colored 
them for a further short burst of computation. Using this "evolve-restrict-rescale-lift" 
procedure repeatedly, the correct functional form of the self-similar solution will arise 
asymptotically: the shape of the coarse-grained description (the ICDF) converges to 
the inverse of the integral of a Gaussian (see Figure 4.) 

In the second numerical experiment, we used a Newton iteration to converge to the 
stationary, dynamically renormalized solution - that is, the self-similar solution. To do 
this, we need to restrict the CDF to a finite (preferably low-dimensional) approximate 
representation. Orthogonal polynomials usually provide computationally simple ap- 
proximations, but unfortunately they are not useful for the CDF which has a possible 
support from — oo to -|-oo (that is, the CDF is c{x) where x has a potentially infinite 
range; c{x) must be monotone non-decreasing and must lie between and 1). This is 
a second, and more important, reason for using the ICDF, x{p) where p lies in [0, 1]. 
This function can easily be approximated by a low-degree polynomial over its range. 

We wish to find a polynomial approximation, Q{p), that approximates the computed 
positions. That is, on the finite set of points, {xi}, corresponding to {pi} = {{i — 
0.5) /N}, we would like Q{pi) ~ Xj. This provides us our restriction of the microscopic 
data. Then we can evaluate Q at any point in [0,1] in the lifting process. (Typically, we 
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will evaluate it at each Pi to get an approximate Xj.) We are interested in minimizing 
the error of the approximation only at the arguments pi,i = 1, ■ ■ ■ , N . If we use least 
squares approximation with weights Wi we thus want to minimize 

N 

- QiPi)fwi. (21) 

i=l 

In the experiments we used unit weights, Wi = 1. The best way to do this com- 
putationally is to create a basis for the TV-dimensional vector space, 0s, s = 1, - ■ ■ ,N 
such that the s-th basis vector, (pg has its i-th element defined as the value of an 
(s — l)-st degree polynomial evaluated at pi, that is, (f)si = QsiPi)- To make this basis 
set orthonormal, we require that 

N 

(0s, 0t) ^ IsipMtiPi) = Sst- (22) 

i=l 

The q polynomials thus defined are scalar multiples of the orthogonal polynomials 
defined in the usual way from the L2 norm over the unit interval. 

Since the steady state solution of w is symmetric, its ICDF will be an odd function 
oi p — 0.5. For this reason, we chose orthogonal polynomials on the finite set {—0.5 + 
{i — 0.5) /N}. For the purposes of the Newton iteration we considered a fifth degree 
polynomial representation of the ICDF written as 

Q{p) = ciMp - 0-5) + C3Mp - 0-5) + c^Mp - 0-5). (23) 

The constancy of the number of colored particles provides one template condition for 
rescaling. In this example, the second template condition was applied by requiring that 
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ci be constant. After each burst of microscopic simulation, the restriction to form Eq. 
(|23|) was performed, and then the polynomial was divided by ci to get a scaled polyno- 
mial with ci = 1. The use of a small number of features of the solution ( "observables" ) 
for the Newton iteration with timesteppers is justified when the long-term dynamics 
of an evolution equation are low-dimensional (i.e. lie on an attracting low-dimensional 
manifold parameterized by a small number of "modes" or "observables", like the q 
here) . 

We can perform variance reduction by using multiple copies. In this case, because 
the box is large enough and the problem is translationally invariant, different parts 
of the box can be used for different realizations of evolving CDFs (and each can be 
colored differently so that they can be distinguished, as long as they are far enough 
apart to be uncorrelated). We take an average of several such realizations (typically 
10, or, for stationary point computations, 100). 

5. Results and discussion 

In order to demonstrate how rescaling can accelerate convergence to the self-similar 
shape we design a piecewise linear CDF with fewer particles in the tail part, as shown 
in Fig. 01 We evolve until t = 300 and rescale, but the tail part of the rescaled CDF 
is still significantly away from the Gaussian. In the right lower corner of Figure \^ 
snapshots of the colored particles in the simulation box at t = 300 and after rescaling 
and lifting are shown. It takes five repetitions of the procedure (evolving and rescaling) 
until the rescaled CDF converges visibly to the Gaussian curve. 

To reduce the effect of noise in estimating functional derivatives or the fixed point 
algorithm, 200 copies of the system at equilibrium are let to run further until t = 1000. 
A perturbation as large as 3% of the initial values of the coarse variables (the Cj) is 
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necessary to ensure a meaningful finite difference estimate of the 2x2 Jacobian of the 
coarse self-similar fixed point computation. As shown in Tabled the fixed point values 
of C3, C5 are cs/ci ^ 0.18, cs/ci 0.075, respectively. 

Comparing the fixed point solution shape (reconstructed based only on the com- 
puted polynomial terms) with the Gaussian, the only visual difference is at the tail 
part, where very few particles exist. The actual solution we have found (the result 
of initializing symmetrically with given values for the first three odd Cj and zero for 
the remaining q, evolving for the given time horizon, and rescaling so that the first 
three odd q have the same values) has now acquired components in the remaining 
odd Cj, and is closer to the Gaussian than its truncation. Using molecular dynamics 
constrained on these three macroscopic observables will give us a better estimate of 
the macroscopic fixed point we located (see (38) for further discussion). Different basis 
sets for the approximation can be used if significant information in the tails is not well 
captured this way. 

The process we have described allows us to compute the shape of the self similar 
solution. If we also need the exponents a and (3 in Eq. ^ we have to know a and b 
defined in Eq. ^ and the value of (\imAT-/A)/{limBT-/B) so that we can apply Eqs. 
© and Q. When an equation is known, a and b are obtained by inspection. When 
it is not known in closed form, ideas similar to those used in |^ can be used to test 
the existence of self-similarity and estimate a and b; in particular, this will involve the 
microscopic simulation for short periods to estimate du/dt in Eq. ^ from multiple 
initial conditions that are scaled (in amplitude and space) versions of each other. 

Once a coarse self-similar solution has been converged to, whether through inte- 
gration or Newton-Raphson of the appropriately renormalized coarse dynamics, we 
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can use again short simulations to estimate At/ A and Bi/B. Starting a microscopic 
simulation at time t and evolving for time the relative values of A and i?, that is, 
A{t + d)/A{t) and B{t + d)/B{t), can be obtained from the scaling needed to re-impose 
the template conditions. These lead to approximations of At/A and Bt/B, from which 
a/ (3 can be estimated using Eq. Q. 

6. Conclusions 

We have described a systematic approach for the computation of coarse self-similar 
solutions in situations where the only available model is a microscopic (in this case 
molecular dynamics) simulator. The procedure uses short, appropriately initialized 
bursts of MD simulation to estimate the coarse timestepper of the template-based 
renormalization flow for the process. In the case we studied here, a 2 x 2 numerical 
Jacobian for the Newton iteration was sufficient; more generally, matrix- free tech- 
niques can be combined with this coarse renormalized timestepper to compute stable 
as well as unstable self-similar solutions, compute "finite times to blow-up" in the ap- 
propriate cases as well as estimate the self-similarity exponents. The stability of the 
self-similar solutions in the co-exploding frame can be probed through the same coarse 
renormalized timestepper and matrix-free eigenanalysis techniques. It is also possible 
to combine the equation- free algorithms presented here with the so-called "gaptooth 
scheme" and "patch dynamics" [HlEnilSlEl- The microscopic computations are not 
performed across all of space, but only in relatively small physical domains ("teeth" 
separated by gaps and connected through appropriate boundary conditions). This ap- 
proach exploits smoothness of the macroscopic observables (e.g. particle density) in 
space as well as time, in order to further reduce the microscopic simulation necessary 
to compute coarse self-similar solutions. 
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It is important to notice that the results are vahd over some scahng regime, while for 
model equations they in principle hold over all scales. In the context of non-Newtonian 
fluid mechanics, such techniques could assist the quantitative detection of self-similarity 
in phenomena such as spreading [41 1142 j . More ambitiously, one may envision the use of 
equation-free coarse renormalization to study the self-similar behavior of the evolution 
of spectra in turbulence studies j43j . 

Extensive research explores the onset of dynamic self-similarity as operating pa- 
rameters vary JTS^ as well as the study of asymptotically self-similar solutions. Self- 
similar solutions may explode or collapse in infinite time (as in the case of the diffusion 
equation we studied here) or can have finite-time blow-ups, in which the solution in- 
creasingly accelerates in time. Conversely it is possible to have self-similar solutions 
which progressively decelerate in time; this may be relevant in the description of glassy 
dynamics [1^ , and it is conceivable that extensions of the approach presented here may 
assist in accelerating the computation of self-similar solutions that progressively slow 
down. 
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Figure Captions 

Fig. 1 a) Time evolution of an initial rectangular density profile by the Id diffusion 
equation, Eq. ©. b) Dynamically renormalized evolution of the same shape using Eq. 
H15|) . c) Cumulative distribution function representation of (b). 

Fig. 2 Rescaling the finite time direct simulation commutes with the dynamic 
renormalization flow. 

Fig. 3 Schematic view of the coarse dynamically renormalized timestepper. 

Fig. 4 Coarse evolution of the cumulative distribution function using coarse renor- 
malized timestepping, starting with a piecewise linear CDF. The inset shows (top) a 
snapshot obtained around the center of the computational domain after t = 300, as 
well as (bottom) the result of restricting, rescaling and lifting this snapshot. 
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Table 1: Dynamically renormalized fixed point computation using Newton-Raphson. 



iteration 


cs/ci 


C5/C1 


|AC3/Ci| + |AC5/Ci| 








1 


1.09626 


1 


0.278266 


0.156855 


0.138427 


2 


0.188089 


0.0816683 


0.01203 


3 


0.176717 


0.0745533 


0.00207797 


4 


0.178675 


0.0747852 


0.000377391 
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Figure 1: a) Time evolution of an initial rectangular density profile by the Id diffusion 
equation, Eq. Q. b) Dynamically renormalized evolution of the same shape using Eq. p5|) . 
c) Cumulative distribution function representation of (b). 
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EVOLVE RESCALED 




Figure 2: Rescaling the finite time direct simulation commutes with the dynamic renormal- 
ization flow. 
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Macroscopic Description 
Figure 3: Schematic view of the coarse dynamically renormalized timestepper. 
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Figure 4: Coarse evolution of the cumulative distribution function using coarse renormal- 
ized timestepping, starting with a piecewise linear CDF. The inset shows (top) a snapshot 
obtained around the center of the computational domain after t = 300, as well as (bottom) 
the result of restricting, rescaling and lifting this snapshot. 
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